Intermittency in passive scalar advection 
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A Lagrangian method for the numerical simulation of the Kraichnan passive scalar model is intro- 
duced. The method is based on Monte-Carlo simulations of tracer trajectories, supplemented by a 
point-splitting procedure for coinciding points. Clean scaling behavior for scalar structure functions 
is observed. The scheme is exploited to investigate the dependence of scalar anomalies on the scaling 
exponent £ of the advecting velocity field. The three-dimensional fourth-order structure function is 
specifically considered. 
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The Kraichnan model of passive scalar advection by a 
self-similar Gaussian white-in-time velocity field [|lj is by 
now a paradigm for intermittency in turbulent systems. 
It was shown by R.H. Kraichnan j|] that, for advect- 
ing velocity fields which are white-in-time (J-correlated) , 
the equal-time correlation functions of the scalar field 9 
obey closed equations of motion. In his 1994 paper ||] 
he used this remarkable property and the so-called "lin- 
ear ansatz" to predict the scaling exponents £n of the 
n-th order scalar structure function S n for all space di- 
mensions d > 2 and for all velocity scaling exponents 
< £ < 2. The predicted £ ra 's are anomalous, for exam- 
ple 2(2 — C4 > 0; in other words the scaling exponents 
cannot be predicted by dimensional analysis. Hence, 
the flatness of scalar increments over a distance r grows 
tx r'' 4-2 '' 2 as r — ► 0, a phenomenon referred to as "in- 
termittency" ||. The linear ansatz was revisited via fu- 
sion rules in Ref. and its validity in the limit £ — > 
was questioned in Ref. || . A different approach was de- 
veloped in Refs. |^-^| in which anomalous scaling has 
its roots in the zero modes (solutions of the unforced 
problem) of the closed exact equations satisfied by the 
scalar correlations 01. Their determination for correla- 
tion functions with more than two or three point presents 
a daunting task which has so far received solutions only 
via perturbation theory. Three limits have been identi- 
fied for the Kraichnan model: large d's A, small £'s |t]] 
and £'s close to the Batchelor limit £ = 2 §,|§. The 
first two expansions are regular, while for the third one 
the relevant small parameter should be y/2 — £. (This is 
due to the preservation of the collinear geometry in the 
Batchelor limit, leading to an angular non- uniformity in 
the perturbation analysis HQ.) 

Numerical simulations have up to now been based on 
the direct integration of the passive scalar partial differ- 
ential equation and have been limited to two dimensions 
J[i| , [L2| . Although the predictions of the linear ansatz 
appear compatible with such simulations, it should be 
noted that such calculations are highly delicate. To wit, 
the difficulty of observing for £2 the known asymptotic 
scaling [Q. 



Our aim here is to propose a different numerical strat- 
egy based on the Monte-Carlo simulation of Lagrangian 
trajectories Jis] ]. For structure functions of finite order 
only a finite number of tracer particles is needed. The 
tracer trajectories are easy to simulate and the calcula- 
tion at each time step only involves a small number of 
random vectors, basically, differences of velocities, rather 
than the whole velocity field. Furthermore, working with 
the tracers naturally allows to measure the scaling of the 
structure functions (r) = ((9(r) — 9(0)) 2n ) vs the in- 
tegral scale L of the forcing. Physically, this means that 
the passive scalar variance injection rate (which equals 
its dissipation rate) and the separation r are kept fixed 
while the integral scale L is varied. In an anomaly-free 
theory, e.g. of the Kolmogorov 1941 type, nothing should 
change in inertial-range statistical quantities. Anomalies 
will here be measured directly through the scaling depen- 
dence on L of the structure functions. 

Specifically, let us consider the passive scalar equation 

d t 6{r, t) + v(r, t) ■ V 0(r, t) = n\7 2 9{r, t) + /(r, t). (1) 

For the Kraichnan model jj] , the velocity and the forcing 
are Gaussian independent processes, both homogeneous, 
stationary, isotropic and white-in-time. The velocity is 
self-similar (no infrared cutoff is needed nor assumed in 
our procedure); the correlations of its increments are 
given by : 



(v a (r,t) vp(r,0)) 



v a (r,t) vp(Q,0)) 



S(t). 



(2) 



As for the forcing, (f{r, t) /(0, 0)) = F (r/L) 5(t), where 
F {r/L) is nearly constant for distances r smaller than 
the integral scale L and decreases rapidly for r ^> L. 

When the molecular diffusivity k is simply ignored, 
and 8 is assumed to vanish in the distant past at t = —T. 
eq. ([!]) can be recast as 9(r, t) — J_ T f(r(s),s)ds, with 
the Lagrangian trajectory defined by the stochastic dif- 
ferential equation dr(s) = v(r(s),s) ds and the final 
condition r(f) = r. Using the Wick rule to calculate 



Gaussian averages over the forcing, the scalar correla- 
tions can be expressed as averages of time integrals of 
F over the statistics of Lagrangian trajectories. Fur- 
thermore, zero-mode ideas suggest the universality of the 
scaling exponents with respect to the choice of F. It is 
then convenient to consider the step function F = 1 for 
r < L and F = for r > L. (The fact that its Fourier 
transform is not positive definite is not relevant for the 
sequel as it actually amounts to taking a complex forcing 
function.) The scalar correlations have then very sim- 
ple expressions, e.g. for the second and the fourth-order 
correlations in the stationary state : 

(0(rr)9(r 2 )) = {T& c] (3) 
(6(n) 6{r 2 )6{r 3 )6{r A )) = + + TftT&) c . 

Here, T^ 2 is the (random) amount of time that two parti- 
cles starting at r% and r 2 and moving backwards in time 
spend with their mutual distance |?"i(s) — f*2(s)| < an d 
(*)c denotes the average over the Lagrangian trajectory 
statistics. Expressions similar to (||) are easily derived 
for higher order correlations. Note that we can exchange 
backward and forward motion in time since, according to 
, the statistical properties of v and —v are the same. 
In the limit n — *■ this procedure, which ignores molec- 
ular diffusion, is actually correct as long as all points 
are distinct. However, if we, e.g., put r\ = r 2 we find 
that (0 2 ), given by (||) is incorrect: it diverges oc T as 
T — > oo. With coinciding points, the correct procedure is 
the point-splitting : the tracer particles must be initially 
separated by a small distance O(e) and the value of the 
correlation function for coinciding points is given by the 
limit e — > 0. This is finite for any £ < 2 because, even for 
e — > 0, the particles reach an O(l) separation in a finite 
time, on account of the Holder non-smooth nature of the 
velocity field (see, e.g., Ref. [Q for this important prop- 
erty of what may be termed a "Richardson walk" ) . It is 
then easily checked that (9 2 ) coincides with the known 
value at k = of the analytical solution jj| and that for 
£ = 2 its divergence is logarithmic in e, as it should be in 
the Batchelor regime. 

In our simulations, the point-splitting operation is 
most conveniently implemented by keeping a non- 
vanishing amount of "molecular noise" . By this we un- 
derstand that the different particles, in addition to being 
swept along by the velocity field, are undergoing indepen- 
dent Brownian motions with a small diffusivity n. This 
Brownian diffusion is relevant only for interparticle dis- 
tances smaller than the dissipation scale r\ = O (k 1 ^). 
The corresponding stochastic equations of motion for the 
case of In tracer particles are : 

dvi{s) = v (ri(s), s) ds + \/2~k Wi(s) ds, i = 1, . . . , 2n, 
n{t) = r t (4) 

where ((w^ Js) (w^ (s')) = 5^5^8(8 - «') and r< 
is the position at the (final) time t. It can be checked 



that the 2n-th order scalar correlation functions given 
by (0), with T{ and t interpreted as Eulerian coordi- 
nates, satisfy Kraichnan's closed equations (for details, 
see Refs. @||). 

The Lagrangian method defined by (||) and (0) is nu- 
merically implemented as follows. Because of homogene- 
ity only differences in positions and velocities matter and 
we can work with 2n — I particles for the moments of 
order 2n. The 2n-th order structure function S2 n re_ 
quires n + 1 configurations of such particles. For ex- 
ample, S 2 (r) = 2 [(6 2 ) - (0(r)0(O))] . Since the velocity 
field is white-in-time, equations such as (^) could present 
the well known Ito-Stratonovich ambiguity |f6| which 
is however absent as a consequence of incompressibil- 
ity. The tracer positions are updated using the clas- 
sical Euler-Ito scheme of order 1/2 [^6|. Thus, dur- 
ing the time interval As the Lagrangian positions for 
each configuration of tracers (r<) (i = 1, . . . , 2n — 1 and 
a = 1, . . . ,d) axe shifted by VA7 {{X l ) a + (Y l ) a ). Here, 
(Xi) and (Yi) are two sets of (2n — 1) d Gaussian ran- 
dom variables chosen independently at each time step 
and with the appropriate correlations. For example, us- 
ing r*2 — ri, . . . , r 2n — r i as dynamical variables, we have 
((X 1 ) 1 (JC 2 ) 3 ) = (K(r 2 )- W '(r 1 )) 1 K(r 3 )-«'(r 1 )) 3 ) 
and {(Xi)i (^2)1) = 2k. (The w'-field has no time depen- 
dence and the same spatial correlations as the w-field.) 
Individual realizations are stopped when all the interpar- 
ticle distances become larger than ten times the largest 
integral scale of interest L max . The number of realiza- 
tions needed for the results reported below is from one 
to several millions. 

A severe test for the Lagrangian method is provided 
by the second-order structure function £2 , whose expres- 
sion is known analytically [Q. Its behavior being non- 
anomalous, a flat scaling in L should be observed. 
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FIG. 1. The 3-D second-order structure function S2 vs L 
for £ = 0.6. The separation r — 2.7 x I0~ 2 , the diffusivity 
k — 1.115 x 10~ 2 and the number of realizations is 4.5 x 10 6 . 



The structure function £2 measured by the Lagrangian 
method is shown in Fig. |l| for £ = 0.6 and d = 3 (all struc- 
ture functions are plotted in log-log coordinates). The 
measured slope is 10 -3 and the error on the constant is 
3%. (These figures are typical also for other values of 
£ studied.) Two remarks are in order. First, it follows 
from the analytic solution, that the constant-in-L behav- 
ior holds for all r < L, including in the dissipative region. 
Physically, this corresponds to the fact that, as r moves 
down in the dissipative region, the energy flux becomes 
smaller and smaller, but still remains independent of L. 
Second, the asymptotic scaling for L <C r goes over into 
the scaling for (8 2 ), namely L 2 ~^; the transition to the 
constant-in-L behavior around r = L is very sharp, on 
account of the step function chosen for F. 
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FIG. 3. Same curve as in Fig. 2 for £ = 0.9. The param- 
eters are r = 2.7 x 10~ 2 , K = 4.4 x 10~ 4 and the number of 
realizations is 8 x 10 6 . 
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FIG. 2. The 3-D fourth-order structure function S4 vs L 
for £ = 0.2. The separation r — 2.7 x 10~ 2 , the diffusivity 
K = 0.247 and the number of realizations is 15 x 10 6 . 
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We applied our method to the determination of 
the anomalies for the fourth-order structure function 
S4(r;L) oc r 2< > 2 (L/r) 2 ^ 2 ^^ in three dimensions. The 
results are summarized in the curve of the anomalous 
correction 2(2 — C4 vs £ presented in Fig. |^. The three 
plots of Si vs L in Figs. ||, ||| and || indicate that the 
Lagrangian simulations require more and more compu- 
tational effort when £ is decreased from 2 to 0. This is 
due mainly to the fact that the three correlation func- 
tions appearing in the expression of S4 have dominant 
contributions scaling as L 2 ( 2 ~4) and L 2 ~^, but they are 
both cancelled in the combination giving S4. Making 
the subdominant contribution of S4 to emerge requires 
stronger and stronger cancellations as £ decreases. For 
all the cases reported the scaling is quite clean, as also 
confirmed by the analysis of local scaling exponents (on 
octaves ratios), whose fluctuations give the conservative 
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FIG. 4. Same curve as in Fig. 2 for £ = 1.75. The pa- 
rameters are r — 2.7 x 10~ 2 , k = 10 -9 and the number of 
realizations is 1.5 x 10 . 

error bars in Fig. ^. 

The dot-dashed line in Fig. |B| is the first-order pertur- 
bative prediction (4/5)£, obtained in Ref. 0. The dashed 
line is a fit of the form 07 + 67 s / 2 with 7 = 2 — £ (the 
parameters are a = 0.06 and b — 1.13), showing that 
the data are compatible with an expansion in ^7. Note 
that a term cx y/j is ruled out by the Holder inequality 
C4 < 2^2 = 27 p7[ . It is interesting to remark that the 
crossing of the curve in Fig. |5| with the monotonically 
decreasing (in £) linear ansatz prediction occurs around 



£ = 1. This is the point farthest from the two limits £ = 
and £ = 2 which both have strongly nonlocal dynamics, 
suggesting a possible relation between deviations from 
the linear ansatz and locality of the interactions pol . 
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FIG. 5. The anomaly 2(^2 — C4 f° r the fourth-order structure 
function in the three-dimensional Kraichnan model. 

We finally note that the two main ingredients of the 
method reported here have in fact a wider range of ap- 
plicability than the determination of anomalies for the 
Kraichnan model. First, the Lagrangian tracer method 
appears more flexible than the integration of the passive 
scalar partial differential equation. The latter permits in 
principle measurement of all the observables and some- 
how corresponds to an infinite number of tracer particles. 
Changing their number according to which specific corre- 
lation function is being investigated seems however to be 
more economic and convenient and should also be of in- 
terest for analyzing the advection by more realistic flow. 
Second, considering the scaling behavior vs the integral 
scale L, rather than vs the separation r, could be useful 
in many situations where the injection rate can be con- 
trolled; this includes the simulation of Navier-Stokes flow 
with white-in-time forcing. Such a procedure presents 
the advantage of giving direct access to the scaling expo- 
nent anomalies, which are quantitative measurements of 
intermittency. 
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